GreenAmpt Function

FUNCTION GreenAmpt &
!
(storm, rain,  ksat, theta, thetas, thetar, phy, itheta, dt, cuminf) &
!
RESULT (inf)


IMPLICIT NONE

!Arguments with intent in
INTEGER (KIND = short), INTENT(IN) :: storm !!when = 1 first steo of storm event
REAL (KIND = float), INTENT(IN) :: rain !!rainfall rate (m/s)
REAL (KIND = float), INTENT(IN) :: ksat !!saturated hydraulic conductivity (m/s)
REAL (KIND = float), INTENT(IN) :: theta !!volumetric water content (m3/m3)
REAL (KIND = float), INTENT(IN) :: thetas !!saturated volumetric water content (m3/m3)
REAL (KIND = float), INTENT(IN) :: thetar !!residual volumetric water content (m3/m3)
REAL (KIND = float), INTENT(IN) :: phy !!suction head across the wetting front (m)
INTEGER (KIND = short), INTENT(IN) :: dt !!time step (s)
REAL (KIND = float), INTENT(IN)  :: itheta !!initial soil moisture

!Arguments with intent inout:
REAL (KIND = float), INTENT(INOUT) :: cuminf !!cumulative infiltration (m)


!local declarations:
REAL (KIND = float) :: inf !!infiltration rate (m/s)
REAL (KIND = float)  :: deltrz  !! difference between soil moisture at beginning
                                !! of storm event and residual soil moisture
REAL (KIND = float) :: test  !!test value for the cumulative infiltration (m)
REAL (KIND = float) :: ftest !!test value for the cumulative infiltration (m)
INTEGER (KIND = short) :: i !loop
INTEGER (KIND = short) :: maxiter = 100 !maximum number of iteration for 
                                        ! computing cumulative infiltration
REAL (KIND = float) :: tol = 0.00001 !tolerance to compute cumulative infiltration (m)
REAL (KIND = float) :: deltatfirst !!for ponding time (case 3)
REAL (KIND = float) :: Fponding !!cumulative infiltration at ponding (case 3)
REAL (KIND = float) :: deltatpond !! new deltat calculated for ponding (case 3)



!-------------------------end of declarations----------------------------------

!if this is the first step of storm event:  
! initialize variables
IF ( storm == 1) THEN

  deltrz = thetas - theta
  IF (deltrz <= 0.) THEN 
    deltrz = 0.
  END IF 
  
  !assume all precipitation is infiltrated
  cuminf = rain * dt
  inf = rain
  RETURN
END IF

!If soil is saturated then set infiltration to zero
IF (theta >= thetas) THEN
  inf = 0.
  RETURN
!If surface is unsaturated then calculate infiltration  
ELSE
  
  IF (storm == 2) THEN !storm is running 
      
     !calculate potential infiltration rate at the beginning of the interval
     !from the known value of cumulative infiltration (eq 5.4.1)
      inf = ksat * (phy * (thetas-itheta) / cuminf + 1.)
      
      !compare potential infiltration to rainfall rate
      IF (inf <= rain) THEN !ponding occurs throughout the interval (case 1)
      
         test = Ksat*dt !first guess for cumulative infiltration
         DO i = 1, maxiter
             IF (i == maxiter) THEN
                CALL Catch ('warning', 'Infiltration', 'maximum number of iterations &
                             reached while updating cumulative infiltration &
                             with green-ampt case 1' )
             END IF
             ftest = cuminf + ksat*dt + (phy*(thetas-itheta))* &
                     LOG ((test+(phy*(thetas-itheta)))/(cuminf+(phy*(thetas-itheta))))
 
             IF (ABS(ftest-test) <= tol) THEN 
                cuminf = ftest
                !compute infiltration rate from cumulative infiltration
                inf = ksat * (phy * (thetas-itheta) / cuminf + 1.)
                IF (inf > rain) inf = rain
                EXIT
             ELSE 
               test = ftest
             END IF 
         END DO 
               
      ELSE !no ponding at the beginning of the interval
      
        !calculate tentative values Assuming that all the rain
        !water is penetrated into the soil      
        test = cuminf + rain * dt
        inf = ksat * (phy * (thetas-itheta) / test + 1.)
        
        IF (inf <= rain) THEN !ponding occurs during the interval (case 3)
           Fponding = ksat * phy * (thetas-itheta) / (rain - ksat)
           deltatfirst = (Fponding - cuminf) / rain
           deltatpond = dt - deltatfirst
           
           !compute cumulative infiltration by substituting cuminf = Fponding
           !and dt = deltatpond
           test = Fponding !Ksat*dt !first guess for cumulative infiltration
           DO i = 1, maxiter
               IF (i == maxiter) THEN
                 CALL Catch ('warning', 'Infiltration', 'maximum number of iterations &
                             reached while updating cumulative infiltration &
                              with green-ampt case 3' )
               END IF
                 ftest = Fponding + ksat*deltatpond + (phy*(thetas-itheta))* &
                     LOG ((test+(phy*(thetas-itheta)))/(Fponding+(phy*(thetas-itheta))))
 
               IF (ABS(ftest-test) <= tol) THEN 
                  cuminf = ftest
                  !compute infiltration rate from cumulative infiltration
                  inf = ksat * (phy * (thetas-itheta) / cuminf + 1.)
                  IF (inf > rain) inf = rain
                  EXIT
               ELSE 
                 test = ftest
               END IF 
            END DO
            
               
        ELSE !no ponding throghout the interval (case 2)
           cuminf = test
           inf = rain
        END IF
        
      END IF
  
  ELSE
     inf = 0.
  END IF
END IF  

END FUNCTION GreenAmpt